Comprehensive genome‑wide analysis of the chicken heat shock protein family: identification, genomic organization, and expression profiles in indigenous chicken with highly pathogenic avian influenza infection

Background Heat shock proteins (HSPs) function as molecular chaperones with critical roles in chicken embryogenesis, immune response to infectious diseases, and response to various environmental stresses. However, little is known on HSP genes in chicken. In this study, to understand the roles of chicken HSPs, we performed genome-wide identification, expression, and functional analyses of the HSP family genes in chicken. Results A total of 76 HSP genes were identified in the chicken genome, which were further classified into eight distinct groups (I-VIII) based on phylogenetic tree analysis. The gene-structure analysis revealed that the members of each clade had the same or similar exon-intron structures. Chromosome mapping suggested that HSP genes were widely dispersed across the chicken genome, except in chromosomes 16, 18, 22, 25, 26, and 28–32, which lacked chicken HSP genes. On the other hand, the interactions among chicken HSPs were limited, indicating that the remaining functions of HSPs could be investigated in chicken. Moreover, KEGG pathway analysis showed that the HSP gene family was involved in the regulation of heat stress, apoptotic, intracellular signaling, and immune response pathways. Finally, RNA sequencing data revealed that, of the 76 chicken HSP genes, 46 were differentially expressed at 21 different growth stages in chicken embryos, and 72 were differentially expressed on post-infection day 3 in two indigenous Ri chicken lines infected with highly pathogenic avian influenza. Conclusions This study provides significant insights into the potential functions of HSPs in chicken, including the regulation of apoptosis, heat stress, chaperone activity, intracellular signaling, and immune response to infectious diseases.

Comprehensive genome-wide analysis of the chicken heat shock protein family: identification, genomic organization, and expression profiles in indigenous chicken with highly pathogenic avian influenza infection Background Heat shock proteins (HSPs) or stress proteins are a class of molecular chaperones that play important roles in protein folding and assembly, immune response, intracellular transport/sorting of proteins, control of cell-cycle and signaling, protection against stress, and apoptosis [1].HSPs are constitutively produced by cells or induced in response to high temperature, hypoxia, infectious diseases, toxins, and several other forms of stress or stimuli [2].Several HSP genes are utilized as housekeeping proteins with high expression levels in non-stressed cells [2].HSPs play crucial roles in restoring unfolded proteins to their functional three-dimensional structure under stress conditions [1].Based on their molecular weight, HSPs are classified into major families including HSP105/110, HSP90, HSP70, HSP60, HSP40 (DNAJ), HSP10, and small HSPs [1,[3][4][5][6].
HSP70 proteins, first discovered and characterized in the early 1960s [4], function as the central components of the cellular network of molecular chaperones and folding catalysts [7].The diverse biological function of the HSP70 gene family, which are well characterized in chickens such as the assembly of multimeric protein complexes, facilitating the intracellular folding of proteins, protein transportation across membranes, and the regulatory responses to the heat shock [8][9][10].The expression of HSP70 proteins can be rapidly upregulated in response to environmental stress [11,12]; the levels of HSP70 family proteins are dynamically regulated upon exposure to extrinsic (pathogenic and environmental) or intrinsic stimuli (physiological, replicative, or oncogenic) [13][14][15].
J-proteins, also called HSP40, J-domain, or DNAJproteins, were initially characterized in E. coli as 41-kDa HSPs [5].The DNAJ-protein, the defining feature of J-proteins, is a compact tetrahelical domain of approximately 70 amino acids with a highly conserved and functionally critical histidine, proline, and aspartic acid tripeptide (HPD) motif [16].In the HSP70-HSP40 co-chaperone system, the association between HSP70 proteins and their substrates requires the binding of an ATP to the ATPase domain, which is hydrolyzed to change the conformation of the binding domain.Thus, several HSP70 substrates can specifically bind to its least conserved C-terminal with a high affinity.However, because the ATPase activity of HSP70s is extremely weak, the J-domain of HSP40s is needed for activating the ATPase domain of HSP70s [15][16][17][18].HSP90 proteins (inducible HSP-α and constitutive HSP-β) [7,19] are found in the cytosol, endoplasmic reticulum, chloroplast, and mitochondria to assist the folding, intracellular transport, maintenance, and degradation of proteins as well as participating in cellular signaling and immune responses to pathogenic stimuli by regulating signaling pathways such as PI3K/AKT, VEGF, and NF-kB pathways [3,20,21].Small heat shock proteins (sHSPs; 15-30 kDa) are one of the five types of molecular chaperones that confer tolerance to heat stress.They are found in all domains of life [3,22] and function to prevent irreversible protein aggregation by binding, rather than refolding, misfolded proteins [23].
However, some recently reported chicken HSPs are still not well-defined; their gene structures, physicochemical properties, and phylogenetic relationships are unknown.In addition, the expression patterns of chicken HSP genes during highly pathogenic avian influenza virus (HPAI) infection remain elusive.In this study, we investigated and characterized all potential HSP genes in chicken, revealing the chromosomal locations, gene structures, and physicochemical properties of chicken HSPs.We also examined the expression levels of chicken HSP genes at different stages of embryonic growth and under conditions of HPAIV infection in two indigenous Ri chicken lines in Vietnam.

Identification of potential chicken HSPs genes
A total of 76 HSP genes were identified in chicken genome: 13 genes encoding small HSP family proteins; 44 HSP40 family; one HSP60 family; eight HSP70 family; one HSP47 family; three HSP90 family; three HSP105/110 family; and three HSF family (Table 1).Among these 76 HSP genes, 35 were novel genes that have not been included in databases such as chD-NAJA2-4, chDNAJB1-B5, chDNAJC19-23, chHSP2-8, chHSPA4, and chHSPA14 (Table 1).The full-length coding sequences of these genes were found in both transcriptome and genome databases.The predicted molecular weights varied from 10.95 kDa (chHSPE1) to 520.759 kDa (chDNAJC29).The chHSP10 gene had the shortest conserved domain with 102 amino acids, and the longest domain with 4580 amino acids was found in the genes chHSPE1 and chDNAJC29.
We identified 13 genes of the small HSP (sHSP) gene family, which contained at least one intact alpha-crystallin domain (HSP20) as confirmed by the Pfam and SMART tools.The 44 genes of the HSP40 gene family were subclassified based on domain conservation into types I, II, and III.In total, 11 HSP40 genes were classified as HSP40 type II (lack zincfinger-like motifs); 29 as HSP40 type III due to the presence of a single J domain; and four (chDNAJA1, chDNAJA2, chDNAJA3, and chDNAJA4; 397-506 amino acids) as HSP40 type I, because they possessed the four characteristic canonical domains J, Glycinephenylalanine(G/F)-rich region, two zinc-finger like motifs, and the carboxyl-terminal, first observed in E. coli (Table 1).Although 23 out of 36 chicken chromosomes contained chHSP40 genes, with the highest number of HSP40 genes found on chromosome 1, no HSP40 gene was found on chromosomes 16-18, 22, 24-26, 28-32, or W (Table 1).A total of 8 HSP70 genes along with their chromosomal positions are presented in Table 1.The genes chHSPA13 and chHSPA14 were mapped to chromosome 1; chHSPA4 and chHSPA9 to chromosome 13; and chHSPA4L, chHSPA2, chHSPA5, and chHSPA8 to chromosomes 4, 5, 17, and 24, respectively.Similarly, three genes were assigned to the HSP90 gene family; chHSP90AA1, chHSP90AB1, and chHSP90B1 were mapped to chromosomes 5, 3, and 1, respectively (Table 1).Three genes, chHSPH1, chHSP105, and chHSP108, were assigned to HSP105/108 gene family and mapped to chromosome 1.Three heat shock factor proteins of the HSF gene family, chHSF1, chHSF2, and chHSF3, were assigned to chromosomes 2 (chHSF1) and 3 (chHSF2-3;   of the chicken genome completely lack HSP genes.To analyze the evolutionary relationships among HSPs in chicken, the protein sequences of the 76 HSPs were used to construct a phylogenetic tree, providing a new perspective for the classification of chicken HSPs. Based on the phylogenetic analysis, chicken HSPs were divided into eight clades (I-VIII; Fig. 1).

Gene structure, domain organization, and physicochemical analysis of chicken HSPs
The physicochemical parameters indicated that the pI of small HSPs ranged from 4.49 to 11.01, with most members of this family being acidic except for chHSP10 and chHSPE1, which were found to be basic (Table 1).The pI levels of chicken HSP40 proteins ranged from 4.66 (chDNAJC24) to 11.01 (chDNAJC30; Table 1).Overall, 23 HSP40 proteins were acidic, while 21 were basic.In the HSP70 family, all members were acidic, with pI levels ranging from 5.03 (HSPA5) to 6.04 (HSPA13; Table 1).Similar results were obtained for the members of the HSP90 family, with chH-SP90B1 having the lowest pI (4.75) and chHSP90AA1 the highest (4.92).HSP105/108 and HSF proteins were also acidic, with pI values ranging from 4.6 to 4.96 Fig. 1 Phylogenetic analysis of chicken HSP genes.Full-length amino acid sequences of 76 chicken HSP genes were aligned using the BioEdit V.07 tool, and the phylogenetic tree was constructed in MEGA7 using the neighbor-joining method with 1000 bootstrap replicates (Table 1).The HSP47 protein was basic (pI, 8.07), and the chicken HSP60 was acidic (pI, 5.6; Table 1).We used the chicken genome to analyze the gene structures of chicken HSP genes for a better understanding of the evolutionary conservation of this family.The exon-intron structures of the chicken HSP genes are shown in Fig. 2A.The numbers of introns/ exons vary significantly among HSP genes, potentially indicative of functional diversity.However, the paralogous gene pairs derived from phylogenetic analysis shared a similar gene structure, and the number of introns in these genes ranged from none to 30 (Fig. 2A).Based on the number of introns, genes could be divided into two models: model 1 with no introns and model 2 containing more than one intron.Five of the 76 HSP genes, chHSP2, chHSP3, chHSP9, LOC772158, and chDNAJB8, were assigned to model 1.Model 2 was typical and contained the remaining 71 genes.
In order to be able to respond to various forms of stress in a timely manner, genes must be rapidly activated, which would be assisted by composed of functional units of protein.In this study, we used the SMART tool for the identification and annotation of protein domains and for the analysis of protein domain architectures (Fig. 2B).The proteins of HSP40 family usually had the DNAJ_CXXCXGXG domain, which contained four cysteine-rich repeats of the motif CXXCXGXG and was embedded in the N-terminus of DNAJ domain except in cases of chDNAJC29 and chDNAJC14.On the other hand, some proteins in the HSP40 family included multiple tandem tetratricopeptide repeats domain [24], which is a structural motif present in a wide range of proteins and mediates protein-protein interactions.It can couple with various domains to perform diverse functions as exemplified by chDNAJC3 and chDNAJC7 (Fig. 2B).The genes in the HSP70 family contained a single or double Pfam HSP70 domain, which mediates protein interaction with other signaling proteins in response to stress or disease.Finally, the other genes in the HSP family contained at least one Pfam HSP domain, which plays an important role in the activation of HSPs or induction of HSP gene expression (Fig. 2B).
To understand the functional diversification of the chicken HSP gene family, a conserved motif analysis was performed (Fig. 2C).We searched for 19 putative motifs in each gene.In general, HSP gene of the same group shared similar motifs.Motif 1 was present in most chicken HSP genes (70/76 HSPs), and motif 2 was present in 66 HSP genes.Seven HSPs, chDNAJC23, chDNAJC15, chDNAJC19, chDNAJC22, chHSP1, chHSP11, and chHIKESHI, had only one motif (Fig. 2C).Overall, the structure of the HSP proteins was conserved among the HSP family members.

Potential regulatory relationship between chicken HSPs
Several studies have reported that HSPs in mammals are associated with or induce the expression of genes involved in various pathways including apoptosis, heat/ cold stress, cytokine release, and immune response.Analysis of protein-protein interactions mainly helps in determining co-expressed and co-functioning proteins that cooperatively trigger or inhibit particular cellular functions [1,2].Interactions of HSPs with other HSP family members were examined using the STRING program (Fig. 3A).This analysis identified only 31 chicken HSP genes and could not identify the remaining 45 genes.Our results indicated strong interactions among chicken HSPs and identified more than 110 functional partners of chicken HSP genes.Specifically, the interactions between chHSP90B1 and chHSPA5 exhibited the highest scores (0.999 and 0.998, respectively) and those between chDNAJC3 and chDNAJB11 showed a high score (0.987 and 0.974, respectively).However, the major predicted interactions among HSPs were merely derived from text-mining data and experimental evidence, and lacked gene fusions.Therefore, co-expression and co-functional studies among HSPs should be conducted to provide experimental evidence.
For the cellular pathway analysis, we entered the 76 HSPs genes into the KEGG pathway-mapping database.The obtained pathway results were then integrated (Fig. 3B).Most HSP genes were associated with 26 pathways, and the major pathways were Chaperones and folding catalysts, Protein processing in the endoplasmic reticulum, Nucleotide-binding, ATPbinding, Stress response, Influenza A, and Mitochondrial biogenesis (Fig. 3B).Our results indicated that chicken HSP family genes mostly functioned in the regulation of gene expression and immune responses to pathogens.

Expression profiles of HSP genes at different growth stages of chicken
To gain insights into the temporal and spatial expression patterns of chicken HSP genes during embryonic development, the chicken RNA-seq data GSE86592 [25] were used to analyze the expression profiles of HSP genes at 21 different growth stages of chicken early embryo (Fig. 4).A total of 30 HSP genes were not expressed at any of the 21 growth stages.The results showed that 46 out of 76 HSP genes were expressed at all 21 growth stages of early embryonic growth at different levels.For example, chDNAJC3, chDNAJC12, and chHSP5 were upregulated in oocyte S1-3 stages Fig. 2 (A) Phylogenetic analysis and gene structures of chicken HSP genes.The unrooted neighbor-joining (NJ) tree was generated in MEGA7 with parameter settings indicated in Fig. 1 and based on full-length amino acid sequences of chicken HSPs.The solid green boxes represent exons, black boxes represent genes up/downstream, and black lines represent introns.(B) Phylogenetic analysis and domain architectures of chicken HSP genes.The unrooted neighbor-joining (NJ) tree was generated in MEGA7 with parameter settings indicated in Fig. 1 and based on full-length amino acid sequences of HSPs in chicken.(C) Conserved motifs of the chicken HSP genes were identified by MEME (http://meme-suite.org/).Grey lines represent non-conserved sequences, and each motif is indicated by a colored box numbered at the bottom.The lengths of the motifs in each protein are proportional but downregulated at all other growth stages (Fig. 4).Moreover, some groups of HSP genes were upregulated at some stages but downregulated at others.For example, chDNAJB8, chDNAJB9, chDNAJB13, chD-NAJA2, chHSPD1, chHSPA4L, chHSPA8, chHSPH1, chHSPAB1, and chHSPAA1 genes were upregulated at EGKVIII_S2-4 and EGKX_S5-7 but downregulated at other stages of early embryonic growth (Fig. 4).In addition, chHSF2BP, chDNAJC9, chDNAJC18, chD-NAJC27, chDNAJC6, chDNAJB2, and chDNAJB4 genes were upregulated at zygote S4-6, EGKI_S1-S4, EGKIII_S4, EGKIII_S5, EGKVI_S1, EGKVI_S5, and EGKVI_S6 and downregulated at other stages (Fig. 4).Our results indicated that the HSP genes were significantly expressed at early embryonic growth stages in chickens and also suggested that HSP genes were responsive to environmental conditions (Fig. 4).

Expression profiles of HSP genes in indigenous Vietnamese ri chickens infected with HPAIV
To investigate the potential functions of HSP genes in response to HPAI infection, we analyzed HSP gene expression in the lung tissues of Vietnamese indigenous Ri chickens (both resistant and susceptible lines) infected with highly pathogenic H5N1 virus on post-infection days 1 and 3 by RNA-Seq analysis.We obtained 19.7-22.3 and 21.3-22.1 million sequence reads for the Ri resistant and susceptible chicken lines, respectively.Of the reads, 88.5-91.4% and 92.5-92.6%were successfully mapped to the chicken genome in cases of the Ri resistant and susceptible lines, respectively [26].We obtained a total of 16,555 transcripts in HPAIV-infected resistant lines and 16,554 transcripts in susceptible ones.Similarly, we found a total of 16,555 and 16,554 transcripts in uninfected control resistant and susceptible lines, respectively, on day 1 and day 3 (data not shown).Then, differentially expressed transcripts (P ≤ 0.05; fold-change ≥ 2) were identified in HPAIV-infected resistant and susceptible lines compared with respective controls.

Validation of RNA-Seq results by RT-qPCR
To validate the RNA-Seq results, RT-qPCR was performed to analyze the expression levels of nine HSP genes in the lung tissues of the two indigenous Ri chicken lines infected with HPAIV compared with the respective uninfected controls (Fig. 6).The expression levels measured by RT-qPCR were consistent with the RNA-Seq results (correlation R 2 = 0.8903 and 0.9391 for indigenous Ri resistant chickens infected with HPAIV at day 1 and day 3, respectively; and correlation R 2 = 0.8788 and 0.9356 for susceptible chicken line infected with HPAIV at day 1 and day 3, respectively) as shown in Fig. 7. Generally, the accuracy of the RNA-Seq analysis was confirmed by the RT-qPCR.

Discussion
Unlike previous studies focusing on one or several HSP families, we performed a genome-wide comprehensive survey of chicken HSPs, characterized all potential HSPs across the chicken genome, estimated potential regulatory relationships among them, analyzed transcriptional expression of HSPs at different stages of embryonic growth and in response to highly pathogenic avian influenza virus infection.Our analysis revealed a total of 76 genes (13 sHSP, 44 HSP40, 1 HSP47, 1 HSP70, 8 HSP70, 3 HSP90, 3 HSP105/110, and 3 HSF genes) in chicken.Chicken HSPs confer protection against a variety of cellular stressors [13] and are notably involved in cytoskeletal rearrangements [2] and apoptosis [9].In this study, chicken HSP genes were observed to be dispersed over 23/36 chromosomes (Table 1), and similar results have been reported in humans [27], mice [19], fish [28], and plants [17,29], suggesting the chicken HSP genes were conserved among species.Computational analysis of the physicochemical properties of a family of proteins helps in understanding their functions in vitro.In this study, 52 out of 76 chicken HSP proteins were observed to be acidic and 24 were basic, as indicated by their pI levels.These observations might be indicative of functional differences among chicken HSP proteins, as similar findings might suggest possibly different roles [2].In any case, an in vivo assessment of chicken HSPs is necessary to understand their functional roles before making valid conclusions.
The phylogenetic, gene structure, and motif analysis of homologous sequences offer a wealth of information by identifying conserved residues crucial to the functions or structures of related proteins [30,31].Based on the phylogenetic, gene structure, and motif analysis of chicken HSPs, we found that the most closely related members of a family or subfamily exhibit similar exon/intron structures and intron numbers, and these results were also consistent with the characteristics defined in the above phylogenetic analysis.In the terminal branch of the phylogenetic tree, the numbers of exons/introns were very similar in some of the sister pairs.These findings indicated that some intron loss, along with intron gain events, might have occurred during the structural evolution of the chicken HSP genes.This pattern was also revealed by motif analysis; the type, order, and number of motifs in proteins of the same family/sub-family were similar but differed from those in proteins of other families and subfamilies.
Moreover, conservation of the structural architecture of HSP genes in humans, mice, plants, and fish has demonstrated that the N and C termini, though variable in sequence and length, are essential in preventing the misfolding of proteins; these observations were neatly validated by our findings [19,28,32].HSP40 proteins are classified into three groups based on the presence of specific conserved regions.Type I DNAJ (DNAJA) proteins are characterized by four repeats of the CxxCxGxG-type zinc-finger region, a C-terminal region, and an N-terminal J-domain followed by a glycine/phenylalanine (G/F)-rich region.Type II DNAJ (DNAJB) proteins are very similar to DNAJA proteins, except that they lack the CxxCxGxGtype zinc-finger region.Type III DNAJ (DNAJC) proteins are the most diverse group, as they only carry the J-domain.The proteins that contain a J-like domain but lack the critical HPD tripeptide are classified as type IV J-proteins (DNAJD) [5,16,32] and functionally characterized based on their role as co-chaperones in binding and regulating the activities of HSP70s [19].A total of 44 candidate chicken HSP40 genes were identified and found to be scattered across the chicken genome, and this number was higher than that in human genome with 41 HSP genes reported [28].A large number of genes were identified in the HSP40 family that should be added to its functional mediatory role in stabilizing the interaction between HSP70 and many substrates in different cellular components to meet the functional of cellular [19,28].In addition, 23 (out of 44) chicken HSP40 genes appeared to be acidic based on their pI levels, and 21 exhibited basic properties, suggesting functional differences that may be useful in wet-lab experiments.DNAJ/HSP40 family members contain the J domain, facilitating their binding to HSP70s, although other domains that critical to their functioning have been identified [27].A total of 8 members of HSP70 family were also identified in chicken with the invariant residues mostly found in the nucleotide-binding domain, where HSP70 proteins interact with the J domains of HSP40 proteins.Interestingly, all HSP70 proteins were predicted to be acidic with a minimal variation among their isoelectric points.This could suggest a functional similarity among the chicken HSP70s, further confirming the reports of conserved functional properties of HSP70 proteins across species [8,21,27].
In chicken, the HSP proteins play an important roles in apoptosis, heat stress response, environmental adaptation, immune response, and cell death, especially HSP70 proteins, though by associated or activated the signaling pathway or interact between HSPs family genes with or without signaling pathways [9,13,18,33].Our analysis demonstrated strong interactions among chicken HSPs, especially HSP70, HSP40, and HSP90 family genes/proteins.Recent research show that the HSP90 family genes play important roles in the etiology of several autoimmune diseases and various infections [7,8].The results also showed that the function of HSPs and their roles in molecular signaling and regulation of the immune system in chicken should be evaluated both in vitro and in vivo.Moreover, the remaining HSP proteins that have not been identified and characterized functionally in chicken should be investigated in future studies.Furthermore, the KEGG pathway analysis showed that the chicken HSP gene family is involved in apoptotic pathways, immune response to the pathogen, and signaling pathways [8].The results indicated that the main functions of chicken HSPs were associated with the regulation of apoptosis, heat stress response, adaptation to environmental conditions, and immune response to infection.
Chicken RNA-Seq data from the public databases were further explored to dissect the expression profiles of the HSP genes analyzed in this study.Our analysis indicated that 30 HSP genes were not expressed in chicken early embryos at any of the 21 growth stages, while the remaining 46 genes were expressed at all 21 growth stages of chicken early embryos, with individual HSP genes showing different expression patterns during embryonic growth.Recent reports have indicated that HSP90, HSP70, and small HSP25 proteins play important roles in the development of embryonic chicken lens [16,34].Our results suggested that HSP genes play an important role at different growth stages of chicken early embryo.Moreover, 72/76 HSP genes were expressed in the lung tissues of the two indigenous Ri chicken lines (resistant and susceptible line) infected with HPAIV.The mRNA expression levels of DNAJB11, DNAJC13, DNAJC21, DNAJC5B, HSP47, HSP7, HSP9, and HSPB3 were significantly upregulated in lung tissues of HPAIV-infected chickens on day 3 of HPAIV infection, while those of DNAJB13, DNAJC14, HSP10, and HSPA2 mRNA were significantly downregulated.The expression levels of other HPS genes showed little change after the HPAIV infection.Recent reports indicate that HSP70, small HSP, HSF, HSP90, and HSP60 family proteins play important roles in response to heat stress and environmental conditions [9,13,14,23].Our study is the first to analyze the expression of HSP family genes in response to HPAIV infection in chickens, and the results suggest that HSP genes have essential functions in the development of the chicken embryo.The functions of the HSP genes in adult chickens, however, has been not investigated and should be examined in the near future.Furthermore, the remaining HSP proteins not characterized here should be investigated in future studies, especially to delineate the functions of HSP family genes in protection against infectious diseases in chicken.

Conclusion
In this study, a total of 76 HSP genes in chicken were identified in a genome-wide survey and classified into 8 groups based on the phylogenetic tree analysis.Protein-protein interactions between HSP proteins were low, and their functional roles in chicken have not been characterized and should be investigated in a future study.On the other hand, KEGG analysis indicated that the chicken HSP gene family is involved in the regulation of apoptotic, heat stress, and immune response pathways.Finally, HSP genes were found to be differentially expressed in early chicken embryos and two HPAIV-infected indigenous Ri chicken lines.Our results will contribute to our understanding on the biological evolution of HSPs in chicken, with possible implications in animal breeding, and further sheds light on intercontinental chicken adaptation mechanisms.

Identification of the HSP family members in chicken
In this study, we used the National Center for Biotechnology Information (NCBI) eukaryotic genome annotation resource database to search for the genome-annotated chicken HSP family genes.The NCBI chicken genome annotation release 103 (Gallus gallus, NCBI Annotation Release 103) contained several gene/protein isoforms.As multiple isoforms represent one common gene/protein, we selected only the first isoform and corresponding protein for further analysis.On the other hand, DNAJ protein was used as a keyword to search for J-proteins, and all candidate proteins were identified using the SMART database (http://smart.embl-heidelberg.de/).The retrieved sequences were translated using the online ORF finder tool (http://www.ncbi.nlm.nih.gov/gorf/gorf.html).Further, the predicted ORFs were verified using the BLASTP tool against the NCBI nonredundant protein sequence database.A total of 76 putative HSPs were identified in the chicken genome.

Gene structures, domains, and phylogenetic analysis of chicken HSP genes
The gene structures of chicken HSP genes were analyzed using the Gene Structure Display Server (http:// gsds.cbi.pku.edu.cn/).Protein-conserved motifs were identified using the MEME software as previously described [31] with the following parameters: optimum width, 10-60; the number of repetitions, any; maximum number of motifs, 15.The isoelectric point (pI) was calculated using the Compute pI/Mw tool (https://web.expasy.org/compute_pi/).The domain organizations of the HSP family proteins were analyzed using the SMART (http://smart.embl-heidelberg.de/), protein family (Pfam) (http://pfam.xfam.org/),and NCBI Bath Web CD-Search (https://www.ncbi.nlm.nih.gov/cdd)databases.The phylogenetic tree of HSP genes was constructed based on protein sequence alignment of the 76 HSP genes in chicken using the MEGA7 software with 1000 bootstrap resampling.
The chicken HSP genes were classified into different groups based on the topology of the phylogenetic tree.

Highly pathogenic avian influenza (HPAI) virus
Highly pathogenic avian influenza virus (HPAIV), A/ duck/Vietnam/QB1207/2012 (H5N1), was used in this study.The viral isolate was propagated in 10-dayold embryonated chicken eggs at 37 ℃ for 48 h.The allantoic fluid (AF) of the eggs was then harvested, and aliquots of the AF were stored at − 80 ℃ until use, according to the OIE guideline (Chap.3.3.4)[35].The 50% egg infectious dose (EID50) of the influenza virus was determined as previously described [36].Briefly, 10-fold serial dilutions of the virus were prepared in PBS, and 100 μl of each dilution was inoculated into the chorioallantoic cavities of 10-day-old embryonated chicken eggs.The eggs were incubated at 37 ℃ for 96 h.Five eggs were infected with each virus dilution.Harvested AF was tested for HA activity using 0.5% RBC according to the OIE guideline (Chap.3.3.4)[35].Calculation of the EID50 of virus suspension was performed using Reed and Muench mathematical technique [37].

Indigenous Ri chicken infection with HPAIV
The specific-pathogen-free HPAIV resistant and susceptible indigenous Ri chicken lines (4 weeks of age) were purchased from the Poultry Research Centre of the National Institute of Animal Science, Vietnam.We conducted genotype analyses on Mx and BF2 to identify traits related to resistance and susceptibility as previously described [26,38].For the Mx protein, chickens carrying the 631-allele polymorphism A were selected as resistant, while those with the 631-allele polymorphism G were designated as susceptible.Similarly, among the BF2 haplotypes, B21 was associated with resistance, and B13 indicated susceptibility [26].Consequently, chickens possessing the Mx(A)/ BF2(B21) genotype were designated as HPAIV-resistant Ri chickens, whereas chickens with the Mx(G)/ BF2(B13) genotype were categorized as HPAIV-susceptible Ri chickens.Infection of the two Vietnamese indigenous Ri chicken lines (resistant and susceptible) was conducted and performed according to the Animal Diseases-Diagnosis procedure (Ministry of Agriculture and Rural Development Vietnam, TCVN8400-26:2014).The protocols and guidelines were approved by the the Institutional Animal Care and Use Committee of the National Institute of Veterinary Research, Vietnam.The study was carried out in compliance with the ARRIVE guidelines.A total of 15 Ri resistant or susceptible chicken per group received intranasal inoculation of AF containing 10 4 EID50 of A/chicken/ Vietnam/NA-01/2019 (H5N1) in 200 μL solution.
Uninfected resistant or susceptible Ri chickens served as controls.All Ri chickens were housed in singlelayer cages, with each cage containing five Ri chickens.These cages were fitted with two food dispensers and a water trough.To prevent any potential spore interference, the control and treated groups were kept in separate poultry houses.All Ri chickens had free access to their respective diets and water.The chicken coops were cleaned twice daily, ensuring good air circulation.An insulation lamp maintained the indoor temperature around 25˚C.Following the viral infection, the chickens were monitored for clinical signs of disease, and the samples were collected on days one and three after infection, following the WHO Manual on Animal Influenza Diagnosis and Surveillance [39].

Expression profiles of HSP genes in chicken at different growth stages using publicly available RNA-Seq data
To explore the expression profiles of HSP genes in chicken at different growth stages, including the oocyte, zygote, and intrauterine embryos from Eyalgiladi and Kochav stage I (EGK.I) to EGK.X, the public high-throughput RNA-seq data of chicken available in the Gene Expression Omnibus (accession number: GSE86592) submitted by Hwang et al. [25] were analyzed.A decision tree-based classification analysis was performed based on the class labels.Log2 trimmed mean of M value normalized values were used to calculate a gene expression matrix and consider library size in each sample by edgeR [25].Spearman's correlation coefficients were used to generate a distance matrix to characterize the linear relationship between class labels and gene expression.Hierarchical cluster analysis of the genes was performed using Cluster version 4.49 (http://www.bram.org/serf/Clusters.php)and the Java TreeView tool (http://jtreeview.sourceforge.net/).Cluster map analysis of the HSP genes detected at different growth stages was carried out using Euclidean distance.The p values were calculated using the righttailed Fisher's exact test, and those lower than 0.01 were considered to indicate significant differences.

RNA extraction and quality analysis
Total RNA was extracted from lung tissues using the TRIzol RNA extraction kit (Invitrogen, Carlsbad, CA, USA), purified using the RNeasy Mini Kit (Qiagen, Germantown, MD, USA), and treated with DNase I (Promega, Madison, WI, USA) according to the manufacturer's instructions.The RNA concentration and quality were further determined using an Agilent 2100 bioanalyzer (Agilent Technologies, San Diego CA, USA) and a Tecan F2000 microplate reader (Tecan Group Ltd., Männedorf, Switzerland).Samples with RNA integrity values greater than 7 and high-quality RNA (28 S/18S > 1) were used in the experiments [26].

RNA-sequencing (RNA-seq) analysis
Potentially existing sequencing adapters and raw quality bases in the raw reads were trimmed using Skewer ver 0.2.2 [40].The cleaned high-quality reads were mapped to the reference genome using the software STAR ver 2.5 [41].Since the sequencing libraries were prepared using Illumina's strand-specific library preparation kit, the strand-specific library option, --library-type = fr-first strand, was applied in the mapping process.
To quantify the mapped reads, Cufflinks ver 2.2.1 [42,43] with the strand-specific library option, --librarytype = fr-first strand, and other default options were used.The gene annotation of the reference genome gg6 from the UCSC genome (https://genome.ucsc.edu) in GTF format was used to obtain gene models, and the expression values were calculated in Fragments Per Kilobase of transcript per Million fragments mapped (FPKM).The genes differentially expressed between the two selected biological conditions were analyzed by Cuffdiff software in the Cufflinks package [42,44] with the strand-specific library option, --librarytype = fr-first strand, and other default options.To compare the expression profiles among the samples, the normalized expression values of the selected few hundred differentially expressed genes were clustered by in-house R scripts.

Pathways and interaction analysis
The cellular pathways associated with the chicken HSP genes were analyzed using the Kyoto Encyclopedia of Genes and Genomes (KEGG; https://www.kegg.jp/kegg/) pathway-mapping database against organismspecific parameters (Gallus gallus; gga) [45,46].The interactions of chicken HSP proteins were analyzed by the Search Tool for the Retrieval of Interacting Genes/ Proteins (STRING, Version 10), as previously described [47].The STRING database facilitates the analysis of gene/protein interactions in an organismspecific manner using commonly available sources, including the NCBI PubMed literature database.Moreover, functional annotation in the form of GO was subsequently extracted from the database using Blast2GO v2.7.1 (http://www.geneontology.org/).

Quantitative real-time PCR (RT-qPCR) analysis of HSP transcripts
For cDNA synthesis, up to 3 μg of RNA was treated with 1.0 unit of DNase I and 1.0 μL of 10x reaction buffer (Thermo Fisher Scientific, Waltham, MA, USA) and incubated for 30 min at 37 ℃.Subsequently, 1.0 μL of 50 mM EDTA was added, and the mixture was incubated at 65 ℃ for 10 min to inactivate DNase I. Reverse transcription was performed then reverse transcribed using a Maxima First Strand cDNA Synthesis Kit (Thermo Fisher Scientific), according to the manufacturer's recommendations.Primers were designed using the Lasergene software (DNASTAR Inc., Madison, WI, USA) and synthesized by PhuSa Genomic Co. Ltd. (CanTho, Vietnam; Table 2).Quantitative PCR was performed using a 2x Ampigen SYBR® Green Master Mix (Enzo Life Sciences, Farmingdale, NY, USA) according to the manufacturer's instructions on a QuantStudio 5 Real-Time PCR System (Thermo Fisher Scientific).Gene expression levels were calculated using the 2 −ΔΔCt method and normalized to those of chicken GAPDH [48].All RT-qPCR experiments were performed in triplicates.

Statistical analyses
Statistical analysis was performed using the IBM SPSS software (SPSS 25.0 for Windows; SPSS, Chicago, IL, USA).The results are expressed as the mean ± standard error (SE) of three independent experiments for each group (n = 3) and were compared between groups using Duncan's multi-comparison method.

Fig. 4
Fig. 4 Expression profiles of HSP genes at 21 different growth stages of early chicken embryo, including the oocyte, zygote, and intrauterine embryos from Eyal-giladi and Kochav stage I (EGK.I) to EGK.X.Different colors correspond to log 2 -transformed values.Red and blue indicate higher and lower relative abundance, respectively

Fig. 5
Fig. 5 Expression profiles of HSP genes in two indigenous Ri chicken lines infected with HPAIV.RNA-seq data were used to produce heat maps.Red and green indicate high and low expression levels, respectively.(A) Expression profiles of chicken HSP genes in indigenous Ri resistant chickens infected with HPAIV compared to the controls.(B) Expression profiles of chicken HSP genes in indigenous Ri susceptible chickens infected with HPAIV compared to the controls.(C) Expression profiles of chicken HSP family genes in indigenous Ri resistant chickens infected with HPAIV compared to indigenous Ri susceptible chickens infected with HPAIV

Fig. 6
Fig. 6 Validation of nine differentially expressed chicken HSP genes.The two indigenous Ri chicken lines were infected with 10 4 EID50 of HPAIV.Lung tissues were isolated on post-infection days 1 and 3, and the transcript levels of indicated genes were determined by reverse transcription followed by quantitative PCR (RT-qPCR).Gene expression levels were normalized to those of GAPDH.Experiments were performed in triplicates on samples from five chickens.Significant differences in mRNA expression levels between two chicken lines are indicated as follows: *, p < 0.05; **, p < 0.01; and ***, p < 0.001.Error bars indicate SE of technical replicates (triplicates)

Fig. 7
Fig. 7 Significant correlations between the results of RT-qPCR and RNA-Seq analyses of gene expression in the lung of two indigenous Ri chicken lines infected with 10 4 EID50 of HPAIV

Table 2
Primers used for quantitative PCR (qPCR) validations